!
!  Robust minimazition functions
!
!  Copyright © 2011 F.Hroch (hroch@physics.muni.cz)

!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

module rfun

  implicit none

  ! robust estimator functions
  ! Hogg in Launer, Wilkinson: Robustness in Statistics

  ! precision of real numbers
  integer, parameter, private :: rp = selected_real_kind(15)
  
  ! robust functions
  real(rp),parameter,private :: huber_a = 1.345_rp
  real(rp),parameter,private :: hampel_a=1.7_rp,hampel_b=3.4_rp,hampel_c=8.5_rp
  real(rp),parameter,private :: Andrews_a = 2.1_rp
  real(rp),parameter,private :: Tukey_a = 6.0_rp
  real(rp),parameter,private :: pi = 3.141592653589793115997963468544185161590576171875_rp

contains

  !--------------------------------------------------------------
  ! 
  ! Robust functions
  !

  ! Huber

  function huber(x)

    real(rp) :: huber
    real(rp),intent(in) :: x
 
    if( abs(x) < huber_a )then
       huber = x
    elseif( abs(x) >= huber_a )then
       huber = sign(huber_a,x)
    else
       huber = sign(huber_a,x)
    endif

  end function huber


  function dhuber(x)

    real(rp) :: dhuber
    real(rp),intent(in) :: x
  
    if( abs(x) < huber_a )then
       dhuber = 1.0_rp
    elseif( abs(x) >= huber_a )then
       dhuber = 0.0_rp
    else
       dhuber = 0.0_rp
    endif

  end function dhuber

 
  ! Hampell

  function hampel(x)

    real(rp) :: hampel
    real(rp),intent(in) :: x
 
    if( abs(x) < hampel_a )then
       hampel = x
    elseif( hampel_a <= abs(x) .and. abs(x) < hampel_b )then
       hampel = sign(hampel_a,x)
    elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
       hampel = hampel_a*(x - sign(hampel_c,x))/(hampel_b - hampel_c)
    elseif(  abs(x) >= hampel_c )then
       hampel = 0.0_rp
    else
       hampel = 0.0_rp
    endif

  end function hampel

  function dhampel(x)

    real(rp) :: dhampel
    real(rp),intent(in) :: x
  
    if( abs(x) < hampel_a )then
       dhampel = 1.0_rp
    elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
       dhampel = -sign(hampel_a/(hampel_b - hampel_c),x)
    else
       dhampel = 0.0_rp
    endif

  end function dhampel


  ! Andrews

  function andrews(x)

    real(rp) :: andrews
    real(rp),intent(in) :: x
 
    if( abs(x) < andrews_a*pi )then
       andrews = sin(x/andrews_a)
    else
       Andrews = 0.0_rp
    endif

  end function andrews

  function dandrews(x)

    real(rp) :: dandrews
    real(rp),intent(in) :: x
  
    if( abs(x) < andrews_a*pi )then
       dandrews = cos(x/andrews_a)/andrews_a
    else
       dandrews = 0.0_rp
    endif

  end function dandrews

  
  ! Tukey

  function tukey(x)

    real(rp) :: tukey
    real(rp),intent(in) :: x
 
    if( abs(x) < tukey_a )then
       tukey = x*(1.0_rp - (x/tukey_a)**2)**2
    else
       tukey = 0.0_rp
    endif

  end function tukey

  function dtukey(x)

    real(rp) :: dTukey
    real(rp),intent(in) :: x
  
    if( abs(x) < tukey_a )then
       dtukey = (1.0_rp - (x/tukey_a)**2)*(1.0_rp - 5.0_rp*(x/tukey_a)**2)
    else
       dtukey = 0.0_rp
    endif

  end function dTukey


end module rfun
